#include <armadillo>
#ifndef MAXIT
#define MAXIT 100
#endif
using namespace arma;

#ifndef _THRES
#define _THRES 0.0000001
#endif
class LALogit
{
	public:

	LALogit(void)
	{	
		boost::math::normal_distribution<> *d=new boost::math::normal_distribution<>(0,1); 
		_d=d;
		_conv=0;
	}
	~LALogit()
	{
		delete _d;
	}
	void operator()(mat X,double *Y,mat m, mat s)
	{
		double e=1;
		double thres=_THRES;
		int p=X.n_cols;
		mat theta(p,1);
		theta.fill(0);
		mat Ym=vect2mat<double>(Y,X.n_rows,1);
		theta=inv(X.t()*X)*X.t()*Ym;
		mat b(p,1);
		int i=0;
		while(e>thres & MAXIT>i)
		{
			b=theta-inv(Hessian(X,Y,theta,s))*Gradiant(X,Y,theta,s);
			e=abserr(theta-b);	
			theta=b;
			i++;
		}
		if(i>MAXIT){cout << "reached max iteration!";_conv=1;	}
		_m=theta;
		_S=-inv(Hessian(X,Y,theta,s));
		

	}

	inline double abserr(mat X)
	{
		double sum=0;
		int n=X.n_rows;
		for(int i=0;i<n;i++)
		{
			sum+=abs<double>(X(i,0));
		}
		return sum;
	}
	inline mat Gradiant(mat X,double *Y,mat theta,mat s)
	{
		int p=X.n_cols;
		mat sum(p,1);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			xb=as_scalar(X(i,span::all)*theta);
			double gpxb=1/(1+exp(-xb));//pdf(*_d,xb);
			double pxb=gpxb*(1-gpxb);//cdf(*_d,xb);
				
			if(Y[i]==1)
			{	
				mat temp=pxb*(1/gpxb)*X(i,span::all);
				sum=sum+temp.t();				
			}else{
							
				mat temp=pxb*(1/(1-gpxb))*X(i,span::all);
				sum=sum-temp.t();				
			}
		}
		return sum-as_scalar(arma::sum(theta,0))/s(0,0);
		
	}

inline	mat Hessian(mat X,double *Y,mat theta,mat s)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);

			double gpxb=1/(1+exp(-xb));
			double pxb=gpxb*(1-gpxb);
			double dpxb=pxb-2*pxb*gpxb;
			mat xx=X(i,span::all).t()*X(i,span::all);	
			if(Y[i]==1)
			{	
				mat temp=((dpxb*gpxb-pow(pxb,2))/pow(gpxb,2))*xx;
				sum=sum+temp;				
			}else{
							
				mat temp=((dpxb*(1-gpxb)+pow(pxb,2))/pow(1-gpxb,2))*xx;
				sum=sum-temp;				
			}
		}
		return sum-1/s(0,0);

	}
/*
	mat Gradiant(mat X,double *Y,mat theta,mat s)
	{
		int p=X.n_cols;
		mat sum(p,1);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			xb=dot(X(i,span::all),theta);
			if(xb!=xb)
			{
				cout << "cb3";
			}
			double gpxb=1/(1+exp(-xb));
			double pxb=exp(-xb);
			
			//cout << "foo1 "<<foo1;
			//cout << "foo2 "<<foo2;
			mat temp=gpxb*pxb*X(i,span::all);
			if(Y[i]==1)
			{	
				sum=sum+temp.t();				
			}else{
							
				sum=sum+temp.t()+X(i,span::all).t();				
			}
		}
		sum+=-as_scalar(arma::sum(theta,0))/s(0,0);
		return sum;

	}


inline	mat Hessian(mat X,double *Y,mat theta,mat s)
	{
		int p=X.n_cols;
		mat sum(p,p);
		sum.fill(0);
		int n=X.n_rows;
		double xb=0;
		for(int i=0;i<n;i++)
		{
			
			xb=as_scalar(X(i,span::all)*theta);
			double pxb=exp(-xb);
			mat xx=X(i,span::all).t()*X(i,span::all);	
			sum+=-pxb*xx/pow(1+pxb,2);	

		}
		return -sum-1/s(0,0);

	}*/
	mat Get_Sig(void){return _S;}
	mat Get_Mu(void){return _m;}
	int Get_conv(void){return _conv;}
	private:
	mat _m;
	mat _S;
	int _conv;
	boost::math::normal_distribution<> *_d; 
};

